Pattern dynamics of cohesive granular particles 
under a plane shear 



Satoshi Takada and Hisao Hayakawa 

Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 

606-8502 Japan 

Abstract. We perform three dimensional molecular dynamics simulations of cohesive granular particles under a plane shear. 
From the simulations, we found that the granular temperature of the system abruptly decreases to zero after reaching the 
critical temperature, where the characteristic time t c \ is approximately represented by t c \ °c (£ — £ C r) - ^ with the dissipation 
rate £ , the critical dissipation rate £ cr and the exponent (5 ~ 0.8. We also found that there exist a variety types of clusters 
depending on the initial density and the dissipation rate. 
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INTRODUCTION 



MODEL AND SETUP 



The interactions among macroscopic granular parti- 
cles, such as sands and granules, are characterized by 
a repulsive and a dissipative ones. The energy dissipa- 
tion through inelastic collisions causes destabilization of 
a uniform state of the system. It is known that there ap- 
pear clusters in such systems [1, 2], which may be un- 
derstood by the hydrodynamic equations [3, 4, 5, 6]. 
When a shear is applied to the granular system, there 
appear high-density region, called "shear band" in the 
two dimensional systems [7, 8], and various shapes of 
clusters such as 2D plug, 2D wave or 3D wave for 
three dimensional systems [8, 9]. There exist many pa- 
pers to estimate the transport coefficients by using ki- 
netic theory [10, 11, 12, 13, 14, 15] and to analyze 
the pattern dynamics by using continuum mechanics 
[7, 16, 17, 18, 19, 20, 21, 22, 23, 24]. 

On the other hand, fine powders of submicron order, 
such as tonner particles or interstellar dusts, have attrac- 
tive forces like an electrostatic force [25] . The existence 
of the attractive force causes some new features due to 
competition of the gas-liquid phase transition [26, 27] 
and the dissipative structure [7] . For instance, nucleation 
process near equilibrium is well understood [28], but that 
under a shear has not been well understood yet. In this 
paper, we try to characterize the nonequilibrium pattern 
formation of fine powders under a plane shear based on 
the three dimensional molecular dynamics simulation. 

In our previous paper [29], we have fixed the dissi- 
pation rate and clarified how the patterns depend on the 
density and the shear rate. In this paper, on the other 
hand, we mainly focus on the role of the dissipation rate 
to the pattern dynamics. 



The system we consider consists of monodiperse N(= 
10,000) spheres, whose radius is a and mass is m. The 
system is a cubic and its linear size is L. We choose x- 
axis and j-axis as the direction of the shear velocity and 
that of the velocity gradient, respectively. The interaction 
among particles is assumed to be described by the trun- 
cated Lennard- Jones potential 



u L \ ri j) = 



( r < r c) 

(r > r c ) 



(1) 



with the well depth £, where rtj = \rij\ is the distance 
between z-th and y'-th (1 < ij < N) particles and r c is 
cut-off length (in this paper, we use r c = 3 a). For the 
dissipation force, we use 



F™(r ij: Vij) = -C0(cr- |r v |)(vy rij)r w 



(2) 



with the dissipation rate £, where Vij is the relative 
velocity vector of z-th and j-th particles defined by Vjj = 
Tij, 0(r) is the step function which is one for r > and 
zero for otherwise, and r is a unit vector proportional to 
r. We note that £ is a dissipation parameter related to the 
coefficient of restitution e, for example, e = 0.998 for 
£ = 0.1 and e = 0.983 for £ = 1.0, respectively, at the 
temperature T = 1.4e. Thus we are interested in weakly 
dissipative situations. This weak dissipation is necessary 
to reach a steady state. Finally, the force acting on i-th 
particle is given by 

Fi = - £ V,y u (r y ) + £F-(r, 7 .v, 7 ). (3) 

In general, we cannot ignore the boundary effect 
[30, 31, 32, 33, 34]. In order to suppress the boundary 



(a) (b) 
FIGURE 1. Typical time evolution of the system for the di- 
lute case p = 0.0463 (f = 0.01, = 0.0173), (a) transitional 
regime at f = 65, 500 and (b) steady regime at t* = 92, 000. 
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FIGURE 2. Typical time evolution of the granular tempera- 
ture T g for p = 0.0463, f = 0.01 with several £*. 



effect to the system, we adopt the Lees-Edwards peri- 
odic boundary condition [35, 36] with the aid of SLLOD 
algorithm [36, 37]: 



dt 
dPj 

dt 



— + Yyi*x, 

m 



—± = Fi-yp yi e x 



(4) 
(5) 



Here r { {x u y u zi), p t (Pxi,P y hPzi) are the position 
and the momentum of z-th particle, respectively, y is 
the shear rate and e x is the unit vector in x direction. 
This algorithm is known that the system is relaxed to 
a uniform shear state near equilibrium. We equilibrate 
the system without the shear and the dissipation until 
t = 200 (ma 2 / e) 1 / 2 , and then we apply the shear and the 
dissipation to the system. 
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FIGURE 3. The dependence of the clustering time t *j on the 
dissipation rate — Q r for p = 0.0463, where the solid line 
denotes t c \ oc (£ — £ cr ) -4 / 5 . 



RESULTS 

At first, we consider a dilute case p = 0.0463, where 
p is the average density No 3 /L 3 . The system size is 
L = 60a for this density. We set To = l.Oe as the initial 
temperature, which is slightly higher than the critical 
temperature for the equilibrium Lennard- Jones fluid [26, 
27]. We fix the shear rate y and change the average 
density p and the dissipation rate £. From our simulation, 
we found two steady phases: (i) uniformly sheared phase 
and (ii) clustered phase. In this paper, we focus on the (ii) 
phase. In this phase, the shape of clusters is a spherical 
droplet for low density. A typical time evolution of the 
configuration of the system is drawn in Fig. 1, where 
r* = r/o, f = y(mo 2 /£) l /\ C* = ^(ma 2 /e) l l 2 and 
t* = t(e/mo 2 ) 1 / 2 . There appears one big droplet, and 
all particles in the system are finally absorbed in it. This 
result corresponds to the case for small shear rate in our 
previous paper [29] . 

We try to characterize this nucleation process 
and the droplet growth in terms of a simple phys- 
ical argument based on the granular temperature 
T g = (m/3N)J^ =1 \vi - v\ 2 [38, 39], where v = v(r,f) 
is the average velocity as a function of r. The time 



evolution of T g is plotted in Fig. 2, where T* = T g /e. 
From this figure, we found that there exist a critical 
temperature T CY and the the critical dissipation rate £ 
corresponding to T cr . In this density, T cr is approximately 
equal to 0.9£. Here we note that this temperature is 
slightly lower than the equilibrium value [26, 27]. For 
C < Ccr, the energy gain by the shear and the energy loss 
by the dissipation are balanced, which keeps the system 
uniform. On the other hand, for £ > £ cr , at an early 
stage, the system remains uniform and T g decreases 
gradually from the initial temperature Tq to the critical 
temperature T cr as time goes on. Then, after reaching 
r cr , there appears a big droplet (Fig. 1(a)) and it absorbs 
all particles in the system. After reaching T g = T cr , T g 
abruptly decreases to become zero. The time when T g 
begins to decrease corresponds to the time at which the 
nucleation takes place. 

In order to characterize the nucleation process we in- 
troduce the "clustering time" t c \, which is defined by the 
time when T g begins to drop as in Fig. 2. The dependence 
of t c \ on £ is plotted in Fig. 3, where t* { = t c i(e/mo 2 ) 1 / 2 




FIGURE 4. Typical shapes of clusters: (a) 2D plug, (b) 2D 
plate, (c) 2D inverse plate, (d) 2D inverse plug, and (e) uniform 
shear state for dense cases. The densities are p = 0.156, 0.305, 
0.723, 0.780 and 0.899, respectively. The shear rate is fixed at 
f =0.1. 



and Q Y = £ cr (ma 2 / e) 1 / 2 . This is approximately given by 

fcl = «(C-Ccr)- P (C>Ccr), (6) 

where a and f3 are constants and £ cr depends on the shear 
rate. Here we note that £ cr is a fitting parameter, and we 
evaluate j3 ~ 0.8 and Q v = 0.0043, 0.017 and 0.068 for 
f = 0.005, 0.01 and 0.02, respectively. 

Next, we consider various density cases: p = 0.0463, 
0.156, 0.305, 0.723, 0.780 and 0.899. Here, these values 
correspond to L = 60a, 40a, 32a, 24a, 23.4a and 
22.32a, respectively. We set To = 1.4£ as the initial 
temperature, which is higher than the critical temperature 
for all densities, and set 7* = 0. 1 as the shear rate. Typical 
shapes of clusters for various densities are drawn in 
Fig. 1(b), 4(a)-(d). We found various shapes: a droplet 
(Fig. 1(b)), 2D plug (Fig. 4(a)), 2D plate (Fig. 4(b)), 2D 
inverse plate (Fig. 4(c)) and 2D inverse plug (Fig. 4(d)). 
These shapes for dilute densities are similar to those 
obtained by the previous work for macroscopic granular 
particles [8]. For denser region, the role of particles and 
vacancies are exchanged in the dilute region, that is, the 
system has a particle hole symmetry with respect to the 
average density. For example, 2D plug (Fig. 4(a)) and 2D 
inverse plug (Fig. 4(d)) are symmetric, and this is also 
true for 2D plate (Fig. 4(b)) and 2D inverse plate (Fig. 
4(c)), however, we can not observe any inverse droplet in 




FIGURE 5. Clusters for strong dissipation rates: (a) a cluster, 
whose shape is heart-like for p = 0.0463 (f = 0.1, = 3.0) 
and (b) a cluster vertical to z-axis for p = 0.305 (7* = 0.1, 
= 0.1). In both cases, the dissipation rate satisfies £ ^> £ cr . 



our simulations. In addition, for liquid-solid coexistence 
region, there appear no clusters (Fig. 4(e)). 

For strong dissipation rate £ ^> £ cr , the shapes of 
clusters depend on the initial conditions. Typical shapes 
are plotted in Fig. 5. For instance, the shape of a droplet is 
not spherical but heart-like in the dilute case (Fig. 5(a)), 
but the shape strongly depends on the initial condition. 
Indeed, we observe a spherical droplet similar to Fig. 
1(b) even when we start from the identical shear rate and 
the density with a different initial configuration. We can 
observe a cluster vertical to the z-axis in an intermediate 
density (Fig. 5(b)). 

Time evolution of T g is similar to Fig. 2. The critical 
temperatures are approximately T cr = 1.26e, 1.28e and 
1.16e for p = 0.0463, 0.305 and 0.456, respectively. 
The former two temperatures are slightly higher than the 
equilibrium values [26, 27], while the latter is slightly 
lower than the equilibrium one. 

CONCLUSION 

We studied the relationships between the shear rate 
and the dissipation rate in the system for cohesive gran- 
ular particles under a plane shear. For the dilute case, we 
found that the granular temperature decreases rapidly to 
zero after reaching the critical temperature T cr , at which 
this characteristic time t c \ depends on the dissipation rate 
£ as t c \ = a(£ — £cr) - ^, and the exponent j3 is approxi- 
mately given 0.8 for p = 0.0463. 

Depending on the average density, we observed vari- 
ous types of clusters. We also found that the shapes of 
clusters for dilute cases and those of inverse clusters for 
dense cases are almost symmetric with respect to the av- 
erage density of the system. Some shapes of clusters de- 
pend on the initial configuration for strongly dissipative 
cases. 
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